home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / DBio.more / readseq.c < prev    next >
Text File  |  1996-07-05  |  36KB  |  1,140 lines
  1. /* File: readseq.c
  2.  * main() program for ureadseq.c, ureadseq.h
  3.  *
  4.  * Reads and writes nucleic/protein sequence in various
  5.  * formats. Data files may have multiple sequences.
  6.  *
  7.  * Copyright 1990 by d.g.gilbert
  8.  * biology dept., indiana university, bloomington, in 47405
  9.  * e-mail: gilbertd@bio.indiana.edu
  10.  *
  11.  * This program may be freely copied and used by anyone.
  12.  * Developers are encourged to incorporate parts in their
  13.  * programs, rather than devise their own private sequence
  14.  * format.
  15.  *
  16.  * This should compile and run with any ANSI C compiler.
  17.  * Please advise me of any bugs, additions or corrections.
  18.  *
  19.  */
  20.  
  21. const char *title
  22.     = "readSeq (1Feb93), multi-format molbio sequence reader.\n";
  23.  
  24.  /*  History
  25.   27 Feb 90.  1st release to public.
  26.    4 Mar 90.  + Gary Olsen format
  27.               + case change
  28.               * minor corrections to NBRF,EMBL,others
  29.               * output 1 file per sequence for gcg, unknown
  30.               * define -DNOSTR for c-libraries w/o strstr
  31.               - readseq.p, pascal version, becomes out-of-date
  32.   24 May 90.  + Phylip 3.2 output format (no input)
  33.   20 Jul 90.  + Phylip 3.3 output (no input yet)
  34.               + interactive output re-direction
  35.               + verbose progress info
  36.               * interactive help output
  37.               * dropped line no.s on NBRF output
  38.               * patched in HyperGCG XCMD corrections,
  39.                 - except for seq. documentation handling
  40.               * dropped the IG special nuc codes, as IG has
  41.                 adopted the standard IUB codes (now if only
  42.                 everyone would adopt a standard format !)
  43.   11 Oct 90.  * corrected bug in reading/writing of EMBL format
  44.  
  45.   17 Oct 91.  * corrected bug in reading Olsen format
  46.                 (serious-deletion)
  47.   10 Nov 91.  * corrected bug in reading some GCG format files
  48.                 (serious-last line duplicated)
  49.               + add format name parsing (-fgb, -ffasta, ...)
  50.               + Phylip v3.4 output format (== v3.2, sequential)
  51.               + add checksum output to all forms that have document
  52.               + skip mail headers in seq file
  53.               + add pipe for standard input == seq file (with -p)
  54.               * fold in parts of MacApp Seq object
  55.               * strengthen format detection
  56.               * clarify program structure
  57.               * remove fixed sequence size limit (now dynamic, sizeof memory)
  58.               * check and fold in accumulated bug reports:
  59.               *   Now ANSI-C fopen(..,"w") & check open failure
  60.               *   Define -DFIXTOUPPER for nonANSI C libraries that mess
  61.                   up toupper/tolower
  62.               = No command-line changes; callers of readseq main() should be okay
  63.               - ureadseq.h functions have changed; client programs need to note.
  64.               + added Unix and VMS Make scripts, including validation tests
  65.  
  66.    4 May 92.  + added 32 bit CRC checksum as alternative to GCG 6.5bit checksum
  67.                 (-DBIGCHECKSUM)
  68.     Aug 92    = fixed Olsen format input to handle files w/ more sequences,
  69.                 not to mess up when more than one seq has same identifier,
  70.                 and to convert number masks to symbols.
  71.               = IG format fix to understand ^L
  72.  
  73.   25-30 Dec 92
  74.               * revised command-line & interactive interface.  Suggested form is now
  75.                   readseq infile -format=genbank -output=outfile -item=1,3,4 ...
  76.                 but remains compatible with prior commandlines:
  77.                   readseq infile -f2 -ooutfile -i3 ...
  78.               + added GCG MSF multi sequence file format
  79.               + added PIR/CODATA format
  80.               + added NCBI ASN.1 sequence file format
  81.               + added Pretty, multi sequence pretty output (only)
  82.               + added PAUP multi seq format
  83.               + added degap option
  84.               + added Gary Williams (GWW, G.Williams@CRC.AC.UK) reverse-complement option.
  85.               + added support for reading Phylip formats (interleave & sequential)
  86.               * string fixes, dropped need for compiler flags NOSTR, FIXTOUPPER, NEEDSTRCASECMP
  87.               * changed 32bit checksum to default, -DSMALLCHECKSUM for GCG version
  88.  
  89.    1Feb93
  90.               = revert GenBank output to a fixed left number width which 
  91.                other software depends on.
  92.           = fix for MSF input to handle symbols in names
  93.           = fix bug for possible memory overrun when truncating seqs for
  94.         Phylip or Paup formats (thanks Anthony Persechini)
  95.  
  96.  */
  97.  
  98.  
  99.  
  100. /*
  101.    Readseq has been tested with:
  102.       Macintosh MPW C
  103.       GNU gcc
  104.       SGI cc
  105.       VAX-VMS cc
  106.    Any ANSI C compiler should be able to handle this.
  107.    Old-style C compilers barf all over the source.
  108.  
  109.  
  110. How do I build the readseq program if I have an Ansi C compiler?
  111. #--------------------
  112. # Unix ANSI C
  113. # Use the supplied Makefile this way:
  114. %  make CC=name-of-c-compiler
  115. # OR do this...
  116. % gcc readseq.c ureadseq.c -o readseq
  117.  
  118. #--------------------
  119. $!VAX-VMS cc
  120. $! Use the supplied Make.Com this way:
  121. $  @make
  122. $! OR, do this:
  123. $ cc readseq, ureadseq
  124. $ link readseq, ureadseq, sys$library:vaxcrtl/lib
  125. $ readseq :== $ MyDisk:[myacct]readseq
  126.  
  127. #--------------------
  128. # Macintosh Simple Input/Output Window application
  129. # requires MPW-C and SIOW library (from APDA)
  130. # also uses files macinit.c, macinit.r, readseqSIOW.make
  131. #
  132. Buildprogram readseqSIOW
  133.  
  134. #--------------------
  135. #MPW-C v3 tool
  136. C  ureadseq.c
  137. C  readseq.c
  138. link -w -o readseq -t MPST -c 'MPS ' ∂
  139.    readseq.c.o Ureadseq.c.o ∂
  140.     "{Libraries}"Interface.o ∂
  141.     "{Libraries}"ToolLibs.o ∂
  142.     "{Libraries}"Runtime.o ∂
  143.     "{CLibraries}"StdClib.o
  144. readseq -i1 ig.seq
  145.  
  146. # MPW-C with NCBI tools
  147.  
  148. set NCBI "{Boot}@molbio:ncbi:"; EXPORT NCBI
  149. set NCBILIB1  "{NCBI}"lib:libncbi.o; export NCBILIB1
  150. set NCBILIB2  "{NCBI}"lib:libncbiobj.o; export NCBILIB2
  151. set NCBILIB3  "{NCBI}"lib:libncbicdr.o; export NCBILIB3
  152. set NCBILIB4  "{NCBI}"lib:libvibrant.o; export NCBILIB4
  153.  
  154. C  ureadseq.c
  155. C  -d NCBI -i "{NCBI}"include: ureadasn.c
  156. C  -d NCBI -i "{NCBI}"include: readseq.c
  157. link -w -o readseq -t MPST -c 'MPS ' ∂
  158.    ureadseq.c.o ureadasn.c.o readseq.c.o  ∂
  159.     {NCBILIB4} {NCBILIB2} {NCBILIB1} ∂
  160.     "{Libraries}"Interface.o ∂
  161.     "{Libraries}"ToolLibs.o ∂
  162.     "{Libraries}"Runtime.o ∂
  163.     "{CLibraries}"CSANELib.o ∂
  164.     "{CLibraries}"Math.o ∂
  165.     "{CLibraries}"StdClib.o
  166.  
  167. ===========================================================*/
  168.  
  169.  
  170.  
  171. #include <stdio.h>
  172. #include <string.h>
  173. #include <ctype.h>
  174.  
  175. #include "ureadseq.h"
  176.  
  177. #pragma segment readseq
  178.  
  179.  
  180.  
  181. static char inputfilestore[256], *inputfile = inputfilestore;
  182.  
  183. const char *formats[kMaxFormat+1] = {
  184.     " 1. IG/Stanford",
  185.     " 2. GenBank/GB",
  186.     " 3. NBRF",
  187.     " 4. EMBL",
  188.     " 5. GCG",
  189.     " 6. DNAStrider",
  190.     " 7. Fitch",
  191.     " 8. Pearson/Fasta",
  192.     " 9. Zuker (in-only)",
  193.     "10. Olsen (in-only)",
  194.     "11. Phylip3.2",
  195.     "12. Phylip",
  196.     "13. Plain/Raw",
  197.     "14. PIR/CODATA",
  198.     "15. MSF",
  199.     "16. ASN.1",
  200.     "17. PAUP/NEXUS",
  201.     "18. Pretty (out-only)",
  202.     "" };
  203.  
  204. #define kFormCount  30
  205. #define kMaxFormName 15
  206.  
  207. const  struct formatTable {
  208.   char  *name;
  209.   short num;
  210.   } formname[] = {
  211.     {"ig",  kIG},
  212.     {"stanford", kIG},
  213.     {"genbank", kGenBank},
  214.     {"gb", kGenBank},
  215.     {"nbrf", kNBRF},
  216.     {"embl", kEMBL},
  217.     {"gcg", kGCG},
  218.     {"uwgcg", kGCG},
  219.     {"dnastrider", kStrider},
  220.     {"strider", kStrider},
  221.     {"fitch", kFitch},
  222.     {"pearson", kPearson},
  223.     {"fasta", kPearson},
  224.     {"zuker", kZuker},
  225.     {"olsen", kOlsen},
  226.     {"phylip", kPhylip},
  227.     {"phylip3.2", kPhylip2},
  228.     {"phylip3.3", kPhylip3},
  229.     {"phylip3.4", kPhylip4},
  230.     {"phylip-interleaved", kPhylip4},
  231.     {"phylip-sequential", kPhylip2},
  232.     {"plain", kPlain},
  233.     {"raw", kPlain},
  234.     {"pir", kPIR},
  235.     {"codata", kPIR},
  236.     {"asn.1", kASN1},
  237.     {"msf", kMSF},
  238.     {"paup", kPAUP},
  239.     {"nexus", kPAUP},
  240.     {"pretty", kPretty},
  241.   };
  242.  
  243. const char *kASN1headline = "Bioseq-set ::= {\nseq-set {\n";
  244.  
  245. /* GWW table for getting the complement of a nucleotide (IUB codes) */
  246. /*                     ! "#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[ \]^_`abcdefghijklmnopqrstuvwxyz{|}~ */
  247. const char compl[] = " !\"#$%&'()*+,-./0123456789:;<=>?@TVGHNNCDNNMNKNNYRYSAABWNRN[\\]^_`tvghnncdnnmnknnyrysaabwnrn{|}~";
  248.  
  249.  
  250.  
  251. const char *formatstr( short format)
  252. {
  253.   if (format < 1 || format > kMaxFormat) {
  254.     switch (format) {
  255.       case kASNseqentry :
  256.       case kASNseqset   : return formats[kASN1-1];
  257.       case kPhylipInterleave:
  258.       case kPhylipSequential: return formats[kPhylip-1];
  259.       default: return "(unknown)";
  260.       }
  261.     }
  262.   else return formats[format-1];
  263. }
  264.  
  265. int parseformat( char *name)
  266. {
  267. #define kDupmatch  -2
  268.   int   namelen, maxlen, i, match, matchat;
  269.   char  lname[kMaxFormName+1];
  270.  
  271.   skipwhitespace(name);
  272.   namelen = strlen(name);
  273.   if (namelen == 0)
  274.     return kNoformat;
  275.   else if (isdigit(*name)) {
  276.     i = atol( name);
  277.     if (i < kMinFormat | i > kMaxFormat) return kNoformat;
  278.     else return i;
  279.     }
  280.  
  281.   /* else match character name */
  282.   maxlen = min( kMaxFormName, namelen);
  283.   for (i=0; i<maxlen; i++) lname[i] = to_lower(name[i]);
  284.   lname[maxlen]=0;
  285.   matchat = kNoformat;
  286.  
  287.   for (i=0; i<kFormCount; i++) {
  288.     match = strncmp( lname, formname[i].name, maxlen);
  289.     if (match == 0) {
  290.       if (strlen(formname[i].name) == namelen) return (formname[i].num);
  291.       else if (matchat == kNoformat) matchat = i;
  292.       else matchat = kDupmatch; /* 2 or more partial matches */
  293.       }
  294.     }
  295.   if (matchat == kNoformat || matchat == kDupmatch)
  296.     return kNoformat;
  297.   else
  298.     return formname[matchat].num;
  299. }
  300.  
  301.  
  302.  
  303. static void dumpSeqList(char *list, short format)
  304. {
  305.   long i, l, listlen;
  306.   char s[256];
  307.  
  308.   listlen = strlen(list);
  309.   printf("Sequences in %s  (format is %s)\n", inputfile, formatstr(format));
  310.   for (i=0, l=0; i < listlen; i++) {
  311.     if (list[i] == (char)NEWLINE) {
  312.       s[l] = '\0'; l = 0;
  313.       puts(s);
  314.       }
  315.     else if (l < 255)
  316.       s[l++] = list[i];
  317.     }
  318.   putchar('\n');
  319. }
  320.  
  321.  
  322.  
  323. void usage()
  324. {
  325.   short   i, midi;
  326.  
  327.   fprintf(stderr,title);
  328.   fprintf(stderr,
  329.   "usage: readseq [-options] in.seq > out.seq\n");
  330.   fprintf(stderr," options\n");
  331. /* ? add -d[igits] to allow digits in sequence data, &/or option to specify seq charset !? */
  332.   fprintf(stderr, "    -a[ll]         select All sequences\n");
  333.   fprintf(stderr, "    -c[aselower]   change to lower case\n");
  334.   fprintf(stderr, "    -C[ASEUPPER]   change to UPPER CASE\n");
  335.   fprintf(stderr, "    -degap[=-]     remove gap symbols\n");
  336.   fprintf(stderr, "    -i[tem=2,3,4]  select Item number(s) from several\n");
  337.   fprintf(stderr, "    -l[ist]        List sequences only\n");
  338.   fprintf(stderr, "    -o[utput=]out.seq  redirect Output\n");
  339.   fprintf(stderr, "    -p[ipe]        Pipe (command line, <stdin, >stdout)\n");
  340.   fprintf(stderr, "    -r[everse]     change to Reverse-complement\n");
  341.   fprintf(stderr, "    -v[erbose]     Verbose progress\n");
  342.   fprintf(stderr, "    -f[ormat=]#    Format number for output,  or\n");
  343.   fprintf(stderr, "    -f[ormat=]Name Format name for output:\n");
  344.   midi = (kMaxFormat+1) / 2;
  345.   for (i = kMinFormat-1; i < midi; i++)
  346.    fprintf( stderr, "        %-20s      %-20s\n",
  347.     formats[i], formats[midi+i]);
  348.  
  349.   /* new output format options, esp. for pretty format: */
  350.   fprintf(stderr, "     \n");
  351.   fprintf(stderr, "   Pretty format options: \n");
  352.   fprintf(stderr, "    -wid[th]=#            sequence line width\n");
  353.   fprintf(stderr, "    -tab=#                left indent\n");
  354.   fprintf(stderr, "    -col[space]=#         column space within sequence line on output\n");
  355.   fprintf(stderr, "    -gap[count]           count gap chars in sequence numbers\n");
  356.   fprintf(stderr, "    -nameleft, -nameright[=#]   name on left/right side [=max width]\n");
  357.   fprintf(stderr, "    -nametop              name at top/bottom\n");
  358.   fprintf(stderr, "    -numleft, -numright   seq index on left/right side\n");
  359.   fprintf(stderr, "    -numtop, -numbot      index on top/bottom\n");
  360.   fprintf(stderr, "    -match[=.]            use match base for 2..n species\n");
  361.   fprintf(stderr, "    -inter[line=#]        blank line(s) between sequence blocks\n");
  362.  
  363.   /******  not ready yet
  364.   fprintf(stderr, "    -code=none,rtf,postscript,ps   code syntax\n");
  365.   fprintf(stderr, "    -namefont=, -numfont=, -seqfont=font   font choice\n");
  366.   fprintf(stderr, "       font suggestions include times,courier,helvetica\n");
  367.   fprintf(stderr, "    -namefontsize=, -numfontsize=, -seqfontsize=#\n");
  368.   fprintf(stderr, "       fontsize suggestions include 9,10,12,14\n");
  369.   fprintf(stderr, "    -namefontstyle=, -numfontstyle=, -seqfontstyle= style  fontstyle for names\n");
  370.   fprintf(stderr, "       fontstyle options are plain,italic,bold,bold-italic\n");
  371.   ******/
  372. }
  373.  
  374. void erralert(short err)
  375. {
  376.   switch (err) {
  377.     case 0  :
  378.       break;
  379.     case eFileNotFound: fprintf(stderr, "File not found: %s\n", inputfile);
  380.       break;
  381.     case eFileCreate: fprintf(stderr, "Can't open output file.\n");
  382.       break;
  383.     case eASNerr: fprintf(stderr, "Error in ASN.1 sequence routines.\n");
  384.       break;
  385.     case eNoData: fprintf(stderr, "No data in file.\n");
  386.       break;
  387.     case eItemNotFound: fprintf(stderr, "Specified item not in file.\n");
  388.       break;
  389.     case eUnequalSize:  fprintf(stderr,
  390.       "This format requires equal length sequences.\nSequence truncated or padded to fit.\n");
  391.       break;
  392.     case eUnknownFormat: fprintf(stderr, "Error: this format is unknown to me.\n");
  393.       break;
  394.     case eOneFormat: fprintf(stderr,
  395.       "Warning: This format permits only 1 sequence per file.\n");
  396.       break;
  397.     case eMemFull: fprintf(stderr, "Out of storage memory. Sequence truncated.\n");
  398.       break;
  399.     default: fprintf(stderr, "readSeq error = %d\n", err);
  400.       break;
  401.     }
  402. } /* erralert */
  403.  
  404.  
  405. int chooseFormat( boolean quietly)
  406. {
  407.   char  sform[128];
  408.   int   midi, i, outform;
  409.  
  410.     if (quietly)
  411.       return kPearson;  /* default */
  412.     else {
  413.       midi = (kMaxFormat+1) / 2;
  414.       for (i = kMinFormat-1; i < midi; i++)
  415.         fprintf( stderr, "        %-20s      %-20s\n",
  416.                         formats[i], formats[midi+i]);
  417.       fprintf(stderr,"\nChoose an output format (name or #): \n");
  418.       gets(sform);
  419.       outform = parseformat(sform);
  420.       if (outform == kNoformat) outform = kPearson;
  421.       return outform;
  422.       }
  423. }
  424.  
  425.  
  426.  
  427. /* read paramater(s) */
  428.  
  429. boolean checkopt( boolean casesense, char *sopt, const char *smatch, short minword)
  430. {
  431.   long  lenopt, lenmatch;
  432.   boolean result;
  433.   short minmaxw;
  434.  
  435.   lenopt = strlen(sopt);
  436.   lenmatch= strlen(smatch);
  437.   minmaxw= max(minword, min(lenopt, lenmatch));
  438.  
  439.   if (casesense)
  440.     result= (!strncmp( sopt, smatch, minmaxw));
  441.   else
  442.     result= (!Strncasecmp( sopt, smatch, minmaxw ));
  443.   /* if (result) { */
  444.     /* fprintf(stderr,"true checkopt(opt=%s,match=%s,param=%s)\n", sopt, smatch, *sparam); */
  445.   /*  } */
  446.   return result;
  447. }
  448.  
  449.  
  450. #define   kMaxwhichlist  50
  451.  
  452. /* global for readopt(), main() */
  453. boolean   chooseall = false, quietly = false, gotinputfile = false,
  454.           listonly = false, closeout = false, verbose = false,
  455.           manyout = false, dolower = false, doupper = false, doreverse= false,
  456.           askout  = true, dopipe= false, interleaved = false;
  457. short     nfile = 0, iwhichlist=0, nwhichlist = 0;
  458. short     whichlist[kMaxwhichlist+1];
  459. long      whichSeq = 0, outform = kNoformat;
  460. char      onamestore[128], *oname = onamestore;
  461. FILE      *foo = NULL;
  462.  
  463. void resetGlobals()
  464. /* need this when used from SIOW, as these globals are not reinited automatically
  465. between calls to local main() */
  466. {
  467.   chooseall = false; quietly = false; gotinputfile = false;
  468.   listonly = false; closeout = false; verbose = false;
  469.   manyout = false; dolower = false; doupper = false; doreverse= false;
  470.   askout  = true; dopipe= false; interleaved = false;
  471.   nfile = 0; iwhichlist=0; nwhichlist = 0;
  472.   whichSeq = 0; outform = kNoformat;
  473.   oname = onamestore;
  474.   foo = NULL;
  475.  
  476.   gPrettyInit(gPretty);
  477. }
  478.  
  479.  
  480. #define kOptOkay  1
  481. #define kOptNone  0
  482.  
  483. int readopt( char *sopt)
  484. {
  485.   char    sparamstore[256], *sparam= sparamstore;
  486.   short   n, slen= strlen(sopt);
  487.  
  488.   /* fprintf(stderr,"readopt( %s) == ", sopt); */
  489.  
  490.   if (*sopt == '?') {
  491.     usage();
  492.     return kOptNone;   /*? eOptionBad or kOptNone */
  493.     }
  494.  
  495.   else if (*sopt == '-') {
  496.  
  497.     char *cp= strchr(sopt,'=');
  498.     *sparam= '\0';
  499.     if (cp) {
  500.       strcpy(sparam, cp+1);
  501.       *cp= 0;
  502.       }
  503.  
  504.     if (checkopt( false, sopt, "-help", 2)) {
  505.       usage();
  506.       return kOptNone;
  507.       }
  508.  
  509.     if (checkopt( false, sopt, "-all", 2)) {
  510.       whichSeq= 1; chooseall= true;
  511.       return kOptOkay;
  512.       }
  513.  
  514.     if (checkopt( false, sopt, "-colspace", 4)) { /* test before -c[ase] */
  515.       n= atoi( sparam);
  516.       gPretty.spacer = n;
  517.       return kOptOkay;
  518.       }
  519.  
  520.     if (checkopt( true, sopt, "-caselower", 2)) {
  521.       dolower= true;
  522.       return kOptOkay;
  523.       }
  524.     if (checkopt( true, sopt, "-CASEUPPER", 2)) {
  525.       doupper= true;
  526.       return kOptOkay;
  527.       }
  528.  
  529.     if (checkopt( false, sopt, "-pipe", 2)) {
  530.       dopipe= true; askout= false;
  531.       return kOptOkay;
  532.       }
  533.  
  534.     if (checkopt( false, sopt, "-list", 2)) {
  535.       listonly = true; askout = false;
  536.       return kOptOkay;
  537.       }
  538.  
  539.     if (checkopt( false, sopt, "-reverse", 2)) {
  540.       doreverse = true;
  541.       return kOptOkay;
  542.       }
  543.  
  544.     if (checkopt( false, sopt, "-verbose", 2)) {
  545.       verbose = true;
  546.       return kOptOkay;
  547.       }
  548.  
  549.     if (checkopt( false, sopt, "-match", 5)) {
  550.       gPretty.domatch= true;
  551.       if (*sparam >= ' ') gPretty.matchchar= *sparam;
  552.       return kOptOkay;
  553.       }
  554.     if (checkopt( false, sopt, "-degap", 4)) {
  555.       gPretty.degap= true;
  556.       if (*sparam >= ' ') gPretty.gapchar= *sparam;
  557.       return kOptOkay;
  558.       }
  559.  
  560.     if (checkopt( false, sopt, "-interline", 4)) {
  561.       gPretty.interline= atoi( sparam);
  562.       return kOptOkay;
  563.       }
  564.  
  565.     if (checkopt( false, sopt, "-item", 2)) {
  566.       char  *cp = sparam;
  567.       nwhichlist= 0;
  568.       whichlist[0]= 0;
  569.       if (*cp == 0) cp= sopt+2; /* compatible w/ old way */
  570.       do {
  571.         while (*cp!=0 && !isdigit(*cp)) cp++;
  572.         if (*cp!=0) {
  573.           n = atoi( cp);
  574.           whichlist[nwhichlist++]= n;
  575.           while (*cp!=0 && isdigit(*cp)) cp++;
  576.           }
  577.       } while (*cp!=0 && n>0 && nwhichlist<kMaxwhichlist);
  578.       whichlist[nwhichlist++]= 0; /* 0 == stopsign for loop */
  579.       whichSeq= max(1,whichlist[0]); iwhichlist= 1;
  580.       return kOptOkay;
  581.       }
  582.  
  583.     if (checkopt( false, sopt, "-format", 5)) {/* -format=phylip, -f2, -form=phylip */
  584.       if (*sparam==0) { for (sparam= sopt+2; isalpha(*sparam); sparam++) ; }
  585.       outform = parseformat( sparam);
  586.       return kOptOkay;
  587.       }
  588.     if (checkopt( false, sopt, "-f", 2)) { /* compatible w/ -fphylip prior version */
  589.       if (*sparam==0) sparam= sopt+2;
  590.       outform = parseformat( sparam);
  591.       return kOptOkay;
  592.       }
  593.  
  594.     if (checkopt( false, sopt, "-output", 3)) {/* -output=myseq */
  595.       if (*sparam==0) { for (sparam= sopt+3; isalpha(*sparam); sparam++) ; }
  596.       strcpy( oname, sparam);
  597.       foo = fopen( oname, "w");
  598.       if (!foo) { erralert(eFileCreate); return eFileCreate; }
  599.       closeout = true;
  600.       askout = false;
  601.       return kOptOkay;
  602.       }
  603.     if (checkopt( false, sopt, "-o", 2)) {  /* compatible w/ -omyseq prior version */
  604.       if (*sparam==0) sparam= sopt+2;
  605.       strcpy( oname, sparam);
  606.       foo = fopen( oname, "w");
  607.       if (!foo) { erralert(eFileCreate); return eFileCreate; }
  608.       closeout = true;
  609.       askout = false;
  610.       return kOptOkay;
  611.       }
  612.  
  613.     if (checkopt( false, sopt, "-width", 2)) {
  614.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  615.       n= atoi( sparam);
  616.       if (n>0) gPretty.seqwidth = n;
  617.       return kOptOkay;
  618.       }
  619.  
  620.     if (checkopt( false, sopt, "-tab", 4)) {
  621.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  622.       n= atoi( sparam);
  623.       gPretty.tab = n;
  624.       return kOptOkay;
  625.       }
  626.  
  627.     if (checkopt( false, sopt, "-gapcount", 4)) {
  628.       gPretty.baseonlynum = false;
  629.       /* if (*sparam >= ' ') gPretty.gapchar= *sparam; */
  630.       return kOptOkay;
  631.       }
  632.     if (checkopt( false, sopt, "-nointerleave", 8)) {
  633.       gPretty.noleaves = true;
  634.       return kOptOkay;
  635.       }
  636.  
  637.     if (checkopt( false, sopt, "-nameleft", 7)) {
  638.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  639.       n= atoi( sparam);
  640.       if (n>0 && n<50) gPretty.namewidth =  n;
  641.       gPretty.nameleft= true;
  642.       return kOptOkay;
  643.       }
  644.     if (checkopt( false, sopt, "-nameright", 7)) {
  645.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  646.       n= atoi( sparam);
  647.       if (n>0 && n<50) gPretty.namewidth =  n;
  648.       gPretty.nameright= true;
  649.       return kOptOkay;
  650.       }
  651.     if (checkopt( false, sopt, "-nametop", 6)) {
  652.       gPretty.nametop= true;
  653.       return kOptOkay;
  654.       }
  655.  
  656.     if (checkopt( false, sopt, "-numleft", 6)) {
  657.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  658.       n= atoi( sparam);
  659.       if (n>0 && n<50) gPretty.numwidth =  n;
  660.       gPretty.numleft= true;
  661.       return kOptOkay;
  662.       }
  663.     if (checkopt( false, sopt, "-numright", 6)) {
  664.       if (*sparam==0) { for (sparam= sopt+2; !isdigit(*sparam) && *sparam!=0; sparam++) ; }
  665.       n= atoi( sparam);
  666.       if (n>0 && n<50) gPretty.numwidth =  n;
  667.       gPretty.numright= true;
  668.       return kOptOkay;
  669.       }
  670.  
  671.     if (checkopt( false, sopt, "-numtop", 6)) {
  672.       gPretty.numtop= true;
  673.       return kOptOkay;
  674.       }
  675.     if (checkopt( false, sopt, "-numbottom", 6)) {
  676.       gPretty.numbot= true;
  677.       return kOptOkay;
  678.       }
  679.  
  680.     else {
  681.       usage();
  682.       return eOptionBad;
  683.       }
  684.     }
  685.  
  686.   else {
  687.     strcpy( inputfile, sopt);
  688.     gotinputfile = (*inputfile != 0);
  689.     nfile++;
  690.     return kOptOkay;
  691.     }
  692.  
  693.  /* return kOptNone; -- never here */
  694. }
  695.  
  696.  
  697.  
  698.  
  699. /* this program suffers some as it tries to be a quiet translator pipe
  700.    _and_ a noisy user interactor
  701. */
  702.  
  703. /* return is best for SIOW, okay for others */
  704. #ifdef SIOW
  705. #define Exit(a)   return(a)
  706. siow_main( int argc, char *argv[])
  707.  
  708. #else
  709. #define Exit(a)   exit(a)
  710.  
  711. main( int argc, char *argv[])
  712. #endif
  713. {
  714. boolean   closein = false;
  715. short     ifile, nseq, atseq, format, err = 0, seqtype = kDNA,
  716.           nlines, seqout = 0, phylvers = 2;
  717. long      i, skiplines, seqlen, seqlen0;
  718. unsigned long  checksum= 0, checkall= 0;
  719. char      *seq, *cp, *firstseq = NULL, *seqlist, *progname, tempname[256];
  720. char      seqid[256], *seqidptr = seqid;
  721. char      stempstore[256], *stemp = stempstore;
  722. FILE      *ftmp, *fin, *fout;
  723. long      outindexmax= 0, noutindex= 0, *outindex = NULL;
  724.  
  725. #define exit_main(err) {        \
  726.   if (closeout) fclose(fout);   \
  727.   if (closein) fclose(fin);   \
  728.   if (*tempname!=0) remove(tempname);\
  729.   Exit(err); }
  730.  
  731. #define indexout()  if (interleaved) {\
  732.   if (noutindex>=outindexmax) {\
  733.     outindexmax= noutindex + 20;\
  734.     outindex= (long*) realloc(outindex, sizeof(long)*outindexmax);\
  735.     if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }\
  736.     }\
  737.   outindex[noutindex++]= ftell(fout);\
  738.   }
  739.  
  740.  
  741.   resetGlobals();
  742.   foo = stdout;
  743.   progname = argv[0];
  744.   *oname = 0;
  745.   *tempname = 0;
  746.   /* initialize gPretty ?? -- done in header */
  747.  
  748.   for (i=1; i < argc; i++) {
  749.     err= readopt( argv[i]);
  750.     if (err <= 0) exit_main(err);
  751.     }
  752.  
  753.                             /* pipe input from stdin !? */
  754.   if (dopipe && !gotinputfile) {
  755.     int c;
  756.     tmpnam(tempname);
  757.     inputfile = tempname;
  758.     ftmp = fopen( inputfile, "w");
  759.     if (!ftmp) { erralert(eFileCreate); exit_main(eFileCreate); }
  760.     while ((c = getc(stdin)) != EOF) fputc(c, ftmp);
  761.     fclose(ftmp);
  762.     gotinputfile= true;
  763.     }
  764.  
  765.   quietly = (dopipe || (gotinputfile && (listonly || whichSeq != 0)));
  766.  
  767.   if (verbose || (!quietly && !gotinputfile)) fprintf( stderr, title);
  768.   ifile = 1;
  769.  
  770.                             /* UI: Choose output */
  771.   if (askout && !closeout && !quietly) {
  772.     askout = false;
  773.     fprintf(stderr,"\nName of output file (?=help, defaults to display): \n");
  774.     gets(oname= onamestore);
  775.     skipwhitespace(oname);
  776.     if (*oname == '?') { usage(); exit_main(0); }
  777.     else if (*oname != 0) {
  778.       closeout = true;
  779.       foo = fopen( oname, "w");
  780.       if (!foo) { erralert(eFileCreate); exit_main(eFileCreate); }
  781.       }
  782.     }
  783.  
  784.   fout = foo;
  785.   if (outform == kNoformat) outform = chooseFormat(quietly);
  786.  
  787.                           /* set up formats ... */
  788.   switch (outform) {
  789.     case kPhylip2:
  790.       interleaved= false;
  791.       phylvers = 2;
  792.       outform = kPhylip;
  793.       break;
  794.  
  795.     case kPhylip4:
  796.       interleaved= true;
  797.       phylvers = 4;
  798.       outform = kPhylip;
  799.       break;
  800.  
  801.     case kMSF:
  802.     case kPAUP:
  803.       interleaved= true;
  804.       break;
  805.  
  806.     case kPretty:
  807.       gPretty.isactive= true;
  808.       interleaved= true;
  809.       break;
  810.  
  811.     }
  812.  
  813.   if (gPretty.isactive && gPretty.noleaves) interleaved= false;
  814.   if (interleaved) {
  815.     fout = ftmp = tmpfile();
  816.     outindexmax= 30; noutindex= 0;
  817.     outindex = (long*) malloc(outindexmax*sizeof(long));
  818.     if (outindex==NULL) { err= eMemFull; erralert(err); exit_main(err); }
  819.     }
  820.  
  821.                         /* big loop over all input files */
  822.   do {
  823.                         /* select next input file */
  824.     gotinputfile = (*tempname != 0);
  825.     while ((ifile < argc) && (!gotinputfile)) {
  826.       if (*argv[ifile] != '-') {
  827.         strcpy( inputfile, argv[ifile]);
  828.         gotinputfile = (*inputfile != 0);
  829.         --nfile;
  830.         }
  831.       ifile++;
  832.       }
  833.  
  834.     while (!gotinputfile) {
  835.       fprintf(stderr,"\nName an input sequence or -option: \n");
  836.       inputfile= inputfilestore;
  837.  
  838.       gets(stemp= stempstore);
  839.       if (*stemp==0) goto fini;  /* !! need this to finish work during interactive use */
  840.       stemp= strtok(stempstore, " \n\r\t");
  841.       while (stemp) {
  842.         err= readopt( stemp); /* will read inputfile if it exists */
  843.         if (err<0) exit_main(err);
  844.         stemp= strtok( NULL, " \n\r\t");
  845.         }
  846.       }
  847.               /* thanks to AJB@UK.AC.DARESBURY.DLVH for this PHYLIP3 fix: */
  848.               /* head for end (interleave if needed) */
  849.     if (*inputfile == 0) break;
  850.  
  851.     format = seqFileFormat( inputfile, &skiplines, &err);
  852.  
  853.     if (err == 0)  {
  854. #ifdef NCBI
  855.       if (format == kASNseqentry || format == kASNseqset)
  856.         seqlist = listASNSeqs( inputfile, skiplines, format, &nseq, &err);
  857.       else
  858. #endif
  859.         seqlist = listSeqs( inputfile, skiplines, format, &nseq, &err);
  860.       }
  861.  
  862.     if (err != 0)
  863.       erralert(err);
  864.  
  865.     else if (listonly) {
  866.       dumpSeqList(seqlist,format);
  867.       free( seqlist);
  868.       }
  869.  
  870.     else {
  871.                                 /* choose whichSeq if needed */
  872.       if (nseq == 1 || chooseall || (quietly && whichSeq == 0)) {
  873.         chooseall= true;
  874.         whichSeq = 1;
  875.         quietly = true; /* no loop */
  876.         }
  877.       else if (whichSeq > nseq && quietly) {
  878.         erralert(eItemNotFound);
  879.         err= eItemNotFound;
  880.         }
  881.       else if (whichSeq > nseq || !quietly) {
  882.         dumpSeqList(seqlist, format);
  883.         fprintf(stderr,"\nChoose a sequence (# or All): \n");
  884.         gets(stemp= stempstore);
  885.         skipwhitespace(stemp);
  886.         if (to_lower(*stemp) == 'a') {
  887.           chooseall= true;
  888.           whichSeq = 1;
  889.           quietly = true; /* !? this means we don't ask for another file 
  890.                             as well as no more whichSeqs... */
  891.           }
  892.         else if (isdigit(*stemp)) whichSeq= atol(stemp);
  893.         else whichSeq= 1; /* default */
  894.         }
  895.       free( seqlist);
  896.  
  897.       if (false /*chooseall*/) {  /* this isn't debugged yet...*/
  898.         fin = fopen(inputfile, "r");
  899.         closein= true;
  900.         }
  901.  
  902.       while (whichSeq > 0 && whichSeq <= nseq) {
  903.                                 /* need to open multiple output files ? */
  904.         manyout = ((chooseall || nwhichlist>1) && nseq > 1
  905.                   && (outform == kPlain || outform == kGCG));
  906.         if (manyout) {
  907.           if ( whichSeq == 1 ) erralert(eOneFormat);
  908.           else if (closeout) {
  909.             sprintf( stemp,"%s_%d", oname, whichSeq);
  910.             freopen( stemp, "w", fout);
  911.             fprintf( stderr,"Writing sequence %d to file %s\n", whichSeq, stemp);
  912.             }
  913.           }
  914.  
  915.         if (closein) {
  916.           /* !! this fails... skips most seqs... */
  917.           /* !! in sequential read, must count seqs already read from whichSeq ... */
  918.           /* need major revision of ureadseq before we can do this */
  919.           atseq= whichSeq-1; /* ? just make whichSeq = 1; atSeq = 0; to read next up */
  920.           seqidptr= seqid;
  921.           seq = readSeqFp( whichSeq, fin, skiplines, format,
  922.                           &seqlen, &atseq, &err, seqidptr);
  923.           skiplines= 0;
  924.           }
  925.         else {
  926.           atseq= 0;
  927.           seqidptr= seqid;
  928. #ifdef NCBI
  929.           if (format == kASNseqentry || format == kASNseqset) {
  930.             seqidptr= NULL;
  931.             seq = readASNSeq( whichSeq, inputfile, skiplines, format,
  932.                      &seqlen, &atseq, &err, &seqidptr);
  933.             }
  934.           else
  935. #endif
  936.           seq = readSeq( whichSeq, inputfile, skiplines, format,
  937.                           &seqlen, &atseq, &err, seqidptr);
  938.           }
  939.  
  940.  
  941.         if (gPretty.degap) {
  942.           char *newseq;
  943.           long newlen;
  944.           newseq= compressSeq( gPretty.gapchar, seq, seqlen, &newlen);
  945.           if (newseq) {
  946.             free(seq); seq= newseq; seqlen= newlen;
  947.             }
  948.           }
  949.  
  950.         if (outform == kMSF) checksum= GCGchecksum(seq, seqlen, &checkall);
  951.         else if (verbose) checksum= seqchecksum(seq, seqlen, &checkall);
  952.         if (verbose)
  953.           fprintf( stderr, "Sequence %d, length= %d, checksum= %X, format= %s, id= %s\n",
  954.                 whichSeq, seqlen, checksum, formatstr(format), seqidptr);
  955.  
  956.         if (err != 0) erralert(err);
  957.         else {
  958.                                   /* format fixes that writeseq doesn't do */
  959.           switch (outform) {
  960.             case kPIR:
  961.               if (seqout == 0) fprintf( foo,"\\\\\\\n");
  962.               break;
  963.             case kASN1:
  964.               if (seqout == 0) fprintf( foo, kASN1headline);
  965.               break;
  966.  
  967.             case kPhylip:
  968.               if (seqout == 0) {
  969.                 if (!interleaved) {  /*  bug, nseq is for 1st infile only */
  970.                   if (chooseall) i= nseq; else i=1;
  971.                   if (phylvers >= 4) fprintf(foo," %d %d\n", i, seqlen);
  972.                   else fprintf(foo," %d %d YF\n", i, seqlen);
  973.                   }
  974.                 seqlen0 = seqlen;
  975.                 }
  976.               else if (seqlen != seqlen0) {
  977.                 erralert(eUnequalSize);
  978.                 if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0);
  979.                 for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
  980.                 seqlen = seqlen0;
  981.                 seq[seqlen] = 0; 
  982.                 }
  983.               break;
  984.  
  985.             case kPAUP:
  986.               if (seqout == 0) {
  987.                 seqtype= getseqtype(seq, seqlen);
  988.                 seqlen0 = seqlen;
  989.                 }
  990.               else if (seqlen != seqlen0) {
  991.                 erralert(eUnequalSize);
  992.                 if (seqlen < seqlen0) seq = (char *)realloc(seq, seqlen0); 
  993.                 for (i=seqlen; i<seqlen0; i++) seq[i]= gPretty.gapchar;
  994.                 seqlen = seqlen0;
  995.                 seq[seqlen] = 0; 
  996.                 }
  997.               break;
  998.  
  999.             }
  1000.  
  1001.           if (doupper)
  1002.             for (i = 0; i<seqlen; i++) seq[i] = to_upper(seq[i]);
  1003.           else if (dolower)
  1004.             for (i = 0; i<seqlen; i++) seq[i] = to_lower(seq[i]);
  1005.  
  1006.           if (doreverse) {
  1007.             long  j, k;
  1008.             char  ctemp;
  1009.             for (j=0, k=seqlen-1; j <= k; j++, k--) {
  1010.               ctemp = compl[seq[j] - ' '];
  1011.               seq[j] = compl[seq[k] - ' '];
  1012.               seq[k] = ctemp;
  1013.               }
  1014.             }
  1015.  
  1016.           if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq != NULL) {
  1017.             for (i=0; i<seqlen; i++)
  1018.               if (seq[i]==firstseq[i]) seq[i]= gPretty.matchchar;
  1019.             }
  1020.  
  1021.  
  1022.           if (gPretty.isactive && gPretty.numtop && seqout == 0) {
  1023.             gPretty.numline = 1;
  1024.             indexout();
  1025.             (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1026.             gPretty.numline = 2;
  1027.             indexout();
  1028.             (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1029.             gPretty.numline = 0;
  1030.             }
  1031.  
  1032.           indexout();
  1033.           nlines = writeSeq( fout, seq, seqlen, outform, seqidptr);
  1034.           seqout++;
  1035.           }
  1036.  
  1037.         if ((gPretty.isactive || outform==kPAUP) && gPretty.domatch && firstseq == NULL) {
  1038.           firstseq= seq;
  1039.           seq = NULL;
  1040.           }
  1041.         else if (seq!=NULL) { free(seq); seq = NULL; }
  1042.  
  1043. #ifdef NCBI
  1044.        if ( (format == kASNseqentry || format == kASNseqset)
  1045.           && seqidptr && seqidptr!= seqid)
  1046.             free(seqidptr);
  1047. #endif
  1048.         if (chooseall) whichSeq++;
  1049.         else if (iwhichlist<nwhichlist) whichSeq= whichlist[iwhichlist++];
  1050.         else whichSeq= 0;
  1051.         }
  1052.       if (closein) { fclose(fin); closein= false; }
  1053.       }
  1054.     whichSeq  = 0;
  1055.   } while (nfile > 0 || !quietly);
  1056.  
  1057.  
  1058. fini:
  1059.   if (firstseq) { free(firstseq); firstseq= NULL; }
  1060.   if (err || listonly) exit_main(err);
  1061.  
  1062.   if (gPretty.isactive && gPretty.numbot) {
  1063.     gPretty.numline = 2;
  1064.     indexout();
  1065.     (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1066.     gPretty.numline = 1;
  1067.     indexout();
  1068.     (void) writeSeq( fout, seq, seqlen, outform, seqidptr);
  1069.     gPretty.numline = 0;
  1070.     }
  1071.  
  1072.   if (outform == kMSF) {
  1073.     if (*oname) cp= oname; else cp= inputfile;
  1074.     fprintf(foo,"\n %s  MSF: %d  Type: N  January 01, 1776  12:00  Check: %d ..\n\n",
  1075.                   cp, seqlen, checkall);
  1076.     }
  1077.  
  1078.   if (outform == kPAUP) {
  1079.     fprintf(foo,"#NEXUS\n");
  1080.     if (*oname) cp= oname; else cp= inputfile;
  1081.     fprintf(foo,"[%s -- data title]\n\n", cp);
  1082.     /* ! now have header lines for each sequence... put them before "begin data;... */
  1083.     }
  1084.  
  1085.   if (outform==kPhylip && interleaved) {
  1086.     if (phylvers >= 4) fprintf(foo," %d %d\n", seqout, seqlen);
  1087.     else fprintf(foo," %d %d YF\n", seqout, seqlen);
  1088.     }
  1089.  
  1090.   if (interleaved) {
  1091.     /* interleave species lines in true output */
  1092.     /* nlines is # lines / sequence */
  1093.     short iline, j, leaf, iseq;
  1094.     char  *s = stempstore;
  1095.  
  1096.     indexout();  noutindex--; /* mark eof */
  1097.  
  1098.     for (leaf=0; leaf<nlines; leaf++) {
  1099.       if (outform == kMSF && leaf == 1) {
  1100.         fputs("//\n\n", foo);
  1101.         }
  1102.       if (outform == kPAUP && leaf==1) {
  1103.         switch (seqtype) {
  1104.           case kDNA     : cp= "dna"; break;
  1105.           case kRNA     : cp= "rna"; break;
  1106.           case kNucleic : cp= "dna"; break;
  1107.           case kAmino   : cp= "protein"; break;
  1108.           case kOtherSeq: cp= "dna"; break;
  1109.           }
  1110.         fprintf(foo,"\nbegin data;\n");
  1111.         fprintf(foo," dimensions ntax=%d nchar=%d;\n", seqout, seqlen);
  1112.         fprintf(foo," format datatype=%s interleave missing=%c", cp, gPretty.gapchar);
  1113.         if (gPretty.domatch) fprintf(foo," matchchar=%c", gPretty.matchchar);
  1114.         fprintf(foo,";\n  matrix\n");
  1115.         }
  1116.  
  1117.       for (iseq=0; iseq<noutindex; iseq++) {
  1118.         fseek(ftmp, outindex[iseq], 0);
  1119.         for (iline=0; iline<=leaf; iline++)
  1120.           if (!fgets(s, 256, ftmp)) *s= 0;
  1121.         if (ftell(ftmp) <= outindex[iseq+1])
  1122.           fputs( s, foo);
  1123.         }
  1124.  
  1125.       for (j=0; j<gPretty.interline; j++)
  1126.         fputs( "\n", foo);  /* some want spacer line */
  1127.       }
  1128.     fclose(ftmp); /* tmp disappears */
  1129.     fout= foo;
  1130.     }
  1131.  
  1132.   if (outform == kASN1)  fprintf( foo, "} }\n");
  1133.   if (outform == kPAUP)  fprintf( foo,";\n  end;\n");
  1134.  
  1135.   if (outindex != NULL) free(outindex);
  1136.   exit_main(0);
  1137. }
  1138.  
  1139.  
  1140.